Add Jones-Wilkins-Lee equation of state with four mixture closures
(isobaric, Kuhl-Khasainov, p-T equilibrium, Rocflu blend) via
jwl_mix_type. New shared m_jwl module, per-fluid eos selector,
GPU-resident parameter arrays, example cases, and inverse-function
test suite.

Without this, the device copy of jwl_mix_type was uninitialized; non-zero
mixture closures (Kuhl, p-T equil, Rocflu) silently executed the wrong
branch on GPU.

…s (5-eq model)

Wire the JWL EOS into the solver hot path for the five-equation model
(model_eqns = 2, Allaire et al. JCP 2002). All changes are guarded by
jwl_idx > 0 and model_eqns == model_eqns_5eq so non-JWL cases are
unaffected.

m_variables_conversion:
- s_convert_conservative_to_primitive_variables: after s_compute_pressure
  gives the stiffened-gas pressure, override it with s_jwl_mix_pressure_er
  using Y = alpha_rho(jwl_idx)/rho and alpha = alpha(jwl_idx). Without
  this, the solver uses the wrong EOS to recover pressure from conserved
  variables, giving thermodynamically inconsistent primitive fields.
- s_convert_primitive_to_conservative_variables: add JWL branch that
  computes E = rho*e_mix(rho,p,Y,alpha) + 0.5*rho|u|^2 + qv via
  s_jwl_mix_energy_pr instead of the stiffened-gas E = gamma*p + pi_inf
  formula, keeping prim->cons consistent with cons->prim.
- s_compute_speed_of_sound: add JWL branch using
  s_jwl_mixture_sound_speed_squared (frozen mass-weighted rule). Accepts
  an optional alpha_rho_j argument for the correct Y_j = alpha_rho_j/rho;
  falls back to the rho0 proxy alpha_j*rho0/rho when not supplied (e.g.
  avg-state calls).
- Export s_jwl_mixture_sound_speed_squared via the module public list.

m_riemann_solver_hll / m_riemann_solver_hllc:
- Add use m_jwl to both modules.
- In every E_L/E_R reconstruction block: when jwl_idx > 0 and 5-eq,
  compute e_L/e_R from pressure via s_jwl_mix_energy_pr and form
  E = rho*e + 0.5*rho|u|^2 instead of the stiffened-gas gamma*p + pi_inf
  formula. This keeps the Riemann wave-speed estimates thermodynamically
  consistent with the JWL EOS for both HLL and HLLC solvers.
- In m_riemann_solver_hll: pass alpha_rho_j to both s_compute_speed_of_sound
  calls so the sound speed uses the actual phasic partial density for Y_j
  rather than the rho0 proxy.

… case to 5-eq

s_compute_speed_of_sound used 'present(alpha_rho_j) .and. alpha_rho_j == alpha_rho_j'
in a single expression. Fortran does not guarantee short-circuit evaluation of .and.,
so gfortran evaluated alpha_rho_j (null pointer) even when it was absent, causing a
SIGSEGV in every HLLC test. Fix: separate into a nested if-block, compute the
rho0-proxy Y_j first, then conditionally override when alpha_rho_j is present and
not NaN.

Also converts examples/1D_jwl_cj_products from the expensive 6-equation pTg relaxation
model (model_eqns=3, relax_model=6) to the 5-equation model (model_eqns=2) with the
isobaric JWL closure, restoring physical CJ detonation parameters (D=6900 m/s,
CR=1.80) that keep the initial pressure (34.5 GPa) above the JWL cold curve (7.1 GPa).

All four JWL tests now pass: 52A7542F, E3460991, DB95E4F7, 25B9827F.